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1 Introduction 



The subject of Physics of fluids and heat transfer have been well established during the last 
century, and has been intensively studied for a wide range of hydrodynamical problems. Many 
current efficient computational fluid dynamics (CFD) software packages offer great flexibility 
in geometry and material properties. However, two-phase flow problems with moving bound- 
aries still present a major challenge to the current state of computational fluid dynamics. 

Although an interface is a three-dimensional region with a thickness of the order of molecular 
diameters it is conveniently modeled as a two dimensional surface (Cosserat surface). Surface 
tension, being a geometrical property of a two dimensional surface, is proportional to the 
curvature of the interface. It is also a material property of the two adjoined materials, thus 
it is an intrinsic property of the surface. In order to include surface tension properly into the 
modeling we need to formulate an additional balance equation for the surface and add it to 
the volume balance equation. 

Interfacial dynamics for Newtonian surface fluids was first described by Scriven in [Scj using 
tensor notation. However, tensor notation is not part of a typical engineering education. Also, 
there is scarce literature on interface balance equations covering surface tension, most use a 
source term instead. 

The main goal of these notes is to give a review of the equations for two phase flow prob- 
lems with an interface between the two phases in a self-contained way, and, in particular, to 
properly include surface tension into the interface balance equations. 

Balance equations at an interface involve geometrical quantities, such as normal and tangen- 
tial vectors on the interface and mean curvature of the interface. We recall them in section [2l 
For the balance equations and the jump conditions we need Reynolds' transport theorems 
for a material body and for a material body with an internal interface. We review them in 
section 13.11 and 13.21 The kinematics of a moving surface are described in section 13. 3i 

In section 14.11 and 14.21 we recall balance equations for single phase problems and ordinary 
jump conditions for an interface, where the adjacent phases are incompressible Newtonian flu- 
ids. Taking surface tension into account increases the complexity. The normal and tangential 
vectors on the interface can be derived by a simple geometric demonstration. However, to de- 
scribe mean curvature (in a physical, intuitive way) use of tensor notation becomes inevitable. 
The jump condition from section |4?2] does not cover surface tension. 

Surface tension is proportional to the curvature of the interface, being a geometrical prop- 
erty of a surface. It is also a material property of the two related materials, thus it is an 
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intrinsic property of the surface. In order to include surface tension properly into a jump 
condition we need to formulate an additional balance equation for the surface and add it to 
the volume balance equation. This is done in section 14.31 In return it involves kinematic 
relations of the interface and the Reynolds transport theorem for surfaces. They are given in 
section 13.31 and 13. 4[ 

Finally, in section [5l the generic balance equations are applied to mass momentum and 
energy and common simplifications and boundary conditions are discussed. 

Preliminaries In this text we use Einstein summation convention, which states a repeated 
index implies a sum over all possible values for that index. When the index takes only the 
values 1,2, Greek letters are used and when the index takes the values 1,2,3, Latin letters 
are used. A vector in a general curvilinear coordinate system can be referred to the standard 
basis {cj} or to the dual basis {e*} and we can write v = v"^ Si oi v = vie^ . The are called 
contravariant components and the Vi are called covariant components of the vector v. The 
dual basis is defined by • e-' =5/- In a Cartesian coordinate system there is no difference 
between standard and dual basis and the contravariant and covariant coordinates of a vector 
are the same. In a curvilinear coordinate system the basis vectors carry a part of the length 
information of a vector. 

The tensor (outer) product is defined by its action 

{w^ v)w = u{v ■ w) or u{v ^ w) = {u ■ v)w . 

Computing the outer products of the base vectors a second order tensor can be represented 
as T = T'-'ej ® ej = T^'-Si ® = tJe"- ® ej = T^e^ (g) The scalar (dot) product between 
two tensors is defined as 

{vi (g) V2) • [wi (g) 1^2) = (f 1 • wi) (f 2 • W2) . 

The gradient of a vector is defined as 
„ 9 

ox* 

and the divergence of a second order tensor is defined as 

Alternatively one could first take the derivative and then compute the dot or tensor product, 
i.e. = (g) Cfc or V • S = ■ e^. Note that both definitions are possible and are used in 
the literature. Our definition is more common in fluid dynamics, the later definition is more 
common in elasticity. 
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2 Geometry of a moving interface 



In this section we summarize the geometrical properties of the moving interface between two 
phases. First, the formulas for normal and tangential vectors, mean curvature and interface 
velocity are given. Then, using an implicit representation of the surface, the geometrical 
properties of two phases at the interface are computed. Two formulas to compute the normal 
vector and two formulas to compute the mean curvature are given. 

The material in this section is mainly based on [Edj . [Arj . [TrTo| and [Erj . They present a 
formulation of classical continuum physics using tensor calculus. For an introduction in tensor 
calculus we refer to |Scha] . |S107| . |S199] and |De] present transport phenomena (including 
interphase transport phenomena) using moderate tensor notation. [Kiij is a differential ge- 
ometry book written more for mathematicians, |0p| has an easier notation but lack some 



derivations. For advanced calculus we refer to Ap and [Kaj . 



2.1 Curvilinear coordinate systems 

A two-dimensional surface can best be an- 
alyzed by covering the surface with a grid, 
see figure [TJ The grid is obtained by the 
curves where and are constant. The 
position of a point on the surface can be 
given intrinsically in terms of the two curvi- 
linear surface coordinates (or parameters) 
and li^, or extrinsically by a position vector 
to the point. This defines a curvilinear co- 
ordinate system which is not orthogonal in 
general. 




Figure 1: Tangent space on a surface 

Tangential and normal vectors If the interface between the two phases is not stationary the 
position vector to a point on the interface is given in cartesian coordinates as 

x{u^,u'^,t) = X^{u^ ,t) El + x'^ {u^ , , t) 62 + X^{u^ ,t) 63 , (1) 

or in index notation as 

a;(n",t) = x^(u",t)ei , (2) 
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where = u^,u'^. 

A Taylor series expansion of x in the surface variables u" up to the linear term yields the 
total derivative 

xiu"" + dn") - xiu"") = dx=^du'' . 

Along n^-curves, = constant (du^ = 0), ^ defines the tangent vector along this curves; 
similarly ^ defines the tangent vector along the -u^ curve. The derivatives with respect to 
the curvilinear coordinates are called covariant derivatives. The covariant derivatives of a 
position vector 

form the base vectors of a local surface coordinate system. In terms of the covariant base 
vectors the surface metric tensor is defined as 

Oa/? = aa-a,p ■ (4) 

The metric tensor is also called the first fundamental form. The local unit normal vector at 
a point normal to the surface is defined by 

ai X a2 

n = r . (5) 

\ai X a2| 

Dual basis Another set of base vectors (a''}, is defined by the surface Kronecker delta, 
aa-a'^ = ^/ , 

they are called dual base vectors or reciprocal or contravariant base vectors respectively. This 
condition orthogonality relation defines a vector that lies in the plane formed by the vectors 
ai, a2, is perpendicular to 02, forms an acute angle with ai. This is also the definition of the 
gradient Vu^ of the surface coordinate, which is perpendicular to the level surface defined by 
u^{x^,x'^,x^,t) = constant. Similarly the orthogonality relation defines the vector a^. Then 
the dual or contravariant base vectors are given by 

a° = Vn« = e, — . (6) 
ox* 

With the contravariant base vectors the covariant surface metric tensor is 

a"^ = a" • a'^ . (7) 



5 



However, often the dual basis is more conveniently calculated by means of the local unit 
normal vector 



a2 X n 



a 



n X ai 



[ai,a2,n\ [ai,a2,n\ 
[ai,a2,n] = [ai x 02] • n is the scalar triple product. 



(8) 



Orthogonal curvilinear coordinate systems If the base vectors of a curvilinear coordinate 
system are mutually orthogonal (ai • 02 = 0), then the covariant and contravariant metric 
tensor reduce simply to 



an 
022 



all 
a22 



For such orthogonal systems the normalized surface vectors are called self reciprocal, in the 



sense that 



ei 



01 a 

ai 



and 



62 



0-2 
a2 



--. It is convenient to introduce unit vectors 



/an V"22 

Surface gradient The identity tensor is defined by 

J = (g) ai + (g) a2 + n n . 

It possesses the property Ix = x for any x. This relation is also called orthogonality relation. 
By subtracting the part related to the normal vector from I the surface identity tensor is 
defined as 



n(sin = a 



an 



(9) 



Similarly the surface gradient is defined by the projection in normal direction subtracted from 
the gradient 



= IgV = (I - n (g) n) V = V - (n (g) n) V . 

By this we get 



a" (g Or 



a 



9x* 



dx^ 



a 



9n" 



(10) 

(11) 
(12) 
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2.2 Mean curvature 



The mean curvature is proportional to the rate of change of the local normal vector with 
respect to the surface coordinates. 



H 



n 



a 



9n° 



n 



(13) 



2.3 Implicit parameterized surface 



To compute the normal and tangential vectors and the mean curvature of the moving interface 
it is necessary to choose a parametrization of the interface. For example, if a vertical tube 
is assumed where a thin film is flowing down along the inner walls (without waves), then 
the problem has rotational symmetry and the interface between the liquid and the gas can 
be parameterized with the surface coordinates = z and = 1}. To assume an implicit 
parameterized surface is more general. 



Normal and tangential vectors of an implicit parametrized surface Every moving surface 
can be locally described by a real-valued function of two variables and time 

z = h{u, V, t) or implicitly F{u, v,z,t) = z — h{u, v,t) = . 

For convenience we write v} = u, = v. By this parametrization the position vector to any 
point on the surface becomes 

a;(n, v,t) = uei + v e2 + h{u, v, t) 63 . 

The tangential vectors are given by the covariant base vectors 



as 



ai 



dx 

du 





dh 
du. 



0-2 



dx 
dv 




1 

dh 
_ dv . 



(14) 



The local unit normal vector is the cross product of the tangent vectors, scaled by its length. 
Alternatively, if the surface is given hy z = h{u,v,t), the unit normal vector can be obtained 
from the gradient of the implicit function F(u, v, z, t) = 0. Expanding F(u, v, z, t) in a Taylor 
series up to the linear term for the variables u,v,z leads to the total derivative 

9F , 9F , 9F , 
aF = — — du + — — dv + — — dz , 
ou OV oz 



7 



which is zero because F = 0. Writing the total derivative as VF • du = where du = 
[du, dv, dz] shows that the gradient VF is perpendicular to the level surface defined by 
F{u, V, z, t) = 0. The local unit normal vector reads then 



n 



VF 

IVFI 



du 
_dh 
dv 

1 



(15) 



Mean curvature of an implicit parametrized surface The mean curvature is given by 
as 



H 



a 



9n 

9n 



dn 

dv 



Alternatively the mean curvature is often computed more conveniently by means of the first 
and second fundamental form as explained below. For the implicit surface parametrization 
the covariant metric tensor, or first fundamental form reads in matrix form 



ai ■ ai 
a2 ■ ai 



ai ■ a2 
a2 ■ a2 



dh dh 
du dv 
dh dh 1 , /dh\2 



on on 1 _|_ / onvz 
dv du ' \dv ' 

The contravariant metric tensor a""^ is defined as the inverse of the covariant metric tensor, 



a"'^ Uyi^ = 6^p. An element of the inverse of a matrix is given by the transpose of the cofactor 
matrix (denoted by a tilde), divided by the determinant of the matrix. This yields for a (2 x 2) 
matrix 



a/3 



1 



det ail ^22 — ai2 a.21 

and for the implicit surface parametrization 
1 



0^22 
-021 



-ai2 

«11 



dh\2 



dh\2 



_dhdh 
du dv 



dh dh 
du dv 
' dh\2 



Both metric tensors are evidently symmetric. They are also positive definite0 
Next, the second fundamental form is defined as 



n 



or alternatively 



9n/3 



(16) 



^The condition that a matrix is positive definite is that aU upper left determinants are positive, 



1 + (t)^ > and a"a^2 ^ (a^)^ = (1 + + {^f) - 



dh dh 
du dv 



dh dh 
du dv 



> 0. 
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The first equation of (I16p yields for our parametrization 



n 



dai 
du 

du 



n 
n 



ov 
da2 
dv 



n 



1 



du'-' 
dvdu 



dudv 
92/1 
dv'' 



f on \'2 

The second fundamental form is also symmetric but not necessarily positive definite. From 
the first and second fundamental form the shape operator or Weingarten map is defined 



L 



'bn 


bi2 


1 


«22 


-ai2 


b2i 


b22 


det Oq,^ 


-021 


ail 



(17) 



L 



1 



It becomes for the implicit surface parametrization 

^1 _j_ f 9^)2\ _ B^fe dh dh d^h /i (dh\2\ 
V ' ^dv' ' dudv du dv dvdu \ ' ^du' ) 

(l _|_ (dh\2\ _ d^h dh dh d^h /-i ( dh\2\ _ 
~^ \dv) J dv^ du dv dv^ \^ ~^ \ du) J 



d^h 
du-' 

d'^h 
dudv 



d^h dh dh, 
du'' du dv 

d^h dh dh 



dudv du dv. 



with ^/o = yl + (|^)^ + (|^)^- L is the product of a symmetric positive definite matrix and 
a symmetric matrix and the eigenvalues of such a product are all real. 

The two eigenvalues ki and K2 of L are called principal curvatures. The mean curvature 
and the GauB curvature K are defined by 



H 



K 



1 



■ trace L 



1 



2 

detL 



2 

det haf 



Ki + K2 , 



Ki K2 . 



det a^is 

So that finally the mean curvature becomes for our parametrization 

^1 + f^)2\ _29^|Zl9Zl+9!|(i + f^)2\ 
V ^dv ■> / dudv du dv ' dv' \ ^du' / 



H 



d^h 

du'^ 



:fu?+i.fvY+^) 



(18) 
(19) 

(20) 



From the two equations for the mean curvature, (jlSp involves only the covariant base vectors 
and the derivatives of the covariant base vectors, whereas involves also the contravariant 
base vectors and their derivatives, so that is often easier to use (|18|) to compute the mean 
curvature, especially in the case of orthogonal coordinate systems. However, the definition of 
the mean curvature with (|13p is more physically intuitive. 

For the computation of the mean curvature with (jlSp only the covariant base vectors and 
the derivatives of the covariant base vectors need to be computed. If (fT3|) is used to compute 
the mean curvature, the contravariant base vectors and their derivatives will also need to 
be computed. Using the shape operator to compute the mean curvature often simplifies the 
computations, especially in the case of orthogonal coordinate systems. 
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3 Kinematics of bulk fluids and of the moving interface 



In this section we present kinematical relations that are necessary to formulate the balance 
equations at an interface between two phases. Experiments show that a fluid interface is in 
fact a three-dimensional region with a thickness on the micro-scale level. Following [Gij , such 
an interface can be regarded as a two-dimensional dividing surface where the effects of the 
interface on the adjoining bulk phases are represented by surface excess mass, momentum and 
energy. 

First we give kinematical relations for the bulk fluids, then we consider a material volume 
with an internal interface. Next we discuss the kinematical relations of the two-dimensional 
moving interface and in the last section we deal with kinematical relations related to the 
curvature of the interface. 

The first part of this section is mainly based on |Ar| . |Scha| and [TrTo] . For the interface 
related sections see also |S107j . |S199j and [Edj . 

3.1 Kinematics of a material volume 

To clarify terminology we first recall some kinematics of bulk fluids and give the Reynolds 
transport theorem for a material volume. 

From the basic assumption of continuum theory, a body consists of infinitely many particles 
without dimension and no space between them and every particle corresponds to a position 
in space. A particle is represented at a given initial time, by a position vector ^, as shown 
in figure [21 The coordinates of ^ are called material coordinates. At another time the same 
particle is represented by another position vector as a function of the initial position of the 
particle and time 



The coordinates of x are called spatial coordinates. 

The initial position of the particle is taken as a reference configuration. Equation (|2ip 
defines the motion of a particle. Assuming continuous motion and that a particle can not 
occupy two places at the same time the relation is a one-to-one mapping and we can also 
write conversely 
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Physical quantities like density, velocity and temperature, which are functions of space and 
time, are called field variables and they are here denoted by ip. A field variable can also be 
given as a function of particle and time. The representation of a field variable as a function 
of space and time is called spatial (or Euler) representation, that is 



if = ip{x, t) 



or 



if = Lp{$,{x,t),t) . 



The representation of a field variable as a function of particle and time is called material (or 
Lagrange) representation, that is 

= fi^, t) or ip = ip{x{^, t),t) . 

Balance equations of mass, momentum and energy are appropriately described in an Eulerian 
framework. 



Material derivative and velocity Field variables are 
derivatives are partial derivatives. Partial derivatives 

stant are denoted by 9. Partial derivatives where ma^ 

denote with an uppercase The partial derivative with re! 



Dip 

15t 



is called material (or convected) derivative and gives the rate oi ychange which an observer 
moving with the particle would see. The material derivative of apfettion-^veetorjsJih^velocity 




so their 
eld con- 
bnstant we 



(22) 



of a given particle 

jjf ' Figure 2: Moving particle 

Balance equations are given in spatial coordinates. To obtain the material derivative of a field 
variable ip{x{$,,t),t) in spatial variables the chain rule has to be applied 

Dip dip dip Dx* 

(23) 



The material derivative is the local rate of change of a given particle at a given position and 
at a given time plus the convective rate of change related to the moving volume. 



Another common notation for the material derivative is a dot on the variable x. 
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Reynolds transport theorem for a material volume For the derivation of the balance equa- 
tions we need the Reynolds transport theorem for a material volume. A mass conserving 
volume is called material volume (or material body) and here denoted by Vq. It is moving 
with time and deforming in general. A quantity Bq continuously defined over a material vol- 
ume Vb is given by Bq = jy^ ipdV. The rate of change of Bo with respect to time is given 

by 

d^o(t) d f 

ip[x, t) aV . 



dt dt 

Vo{x,t) 

In an Eulerian representation VQ{x,t) depends on time, so that integration and differentiation 
can not be interchanged. With the Jacobian J = det the volume element can be 

transformed from spatial coordinates into material coordinates dV = JdVb. The material 
volume element dVb does not depend on time, so that then integration and differentiation can 
be interchanged. The time derivative becomes the material derivative 



Vo{x,t) Vo(i,t) 



Using the material derivative of the Jacobian ^ = J V • u we get 



where we have dropped the integration limits for simplicity. After transforming the volume 
element back into spatial variables and by using the material derivative (|23p of a field variable 
we get for the rate of change with time of Bq 

^I'^dV = [(^+^V-v]dV, (24) 



dij J \'Dt 

Vb 

^ + V-[ipv]^ dV. 

Note that when the derivative of the integral is taken the integration domain has to be 
indicated. By using Gauss theorenjfl, the divergence term in the volume integral can be 
changed into an area integral 

^ J ifdV = J^dV + j) ifv-ndA, (25) 
Vb 



^The Gauss theorem or divergence theorem for a vector / is given as: /V- / dV = f f ndA. 
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where the normal vector is directed outwards on the surface. The velocity v is the velocity 
of mass while moving across the surface. Equations ()24p - ()25p are called Reynolds transport 
theorem. In the form of (I25p the Reynolds transport theorem has a physical meaning: The 
rate of accumulation of a quantity in a material volume can be interpreted as the rate of 
accumulation of the quantity in a volume that equals the material volume at a given time 
plus convective flux (connected to mass) leaving the volume through the surface at that time. 

3.2 Reynolds transport theorem for a material volume with an interface 

We consider a material volume Vq = Vi + 
Vg, where the field variable has the value tpi 
(/ for liquid phase) in the volume V/ and the 
value ifg {g for gas phase) in the volume Vg, 
as shown in figure O An interface between 
two immiscible fluids is called a material in- 
terface, it is formed by the same material 
elements or particles at all times. If phase 
change occurs at an interface between two 
aggregate states of a fluid, as is the case of 
condensation or evaporation, the surface ve- Figure 3: Material volume with interface 
locity u of the interface differs from the ve- 
locity V of the mass, and the interface is called a singular interface. 

The difference between the two values at the surface is denoted by 

(y?! := Lpg - ipi 

and called the jump of ip across the interface. 

The rate of change of Bq with respect to time is the sum of the rate of change of Bi and Bg 
with respect to time 

dBojt) _ dBijt) dBgjt) 
dt dt dt ' 

that is 

The volumes Vi and Vg are not material, so that we need the Reynolds transport theorem in 
a modified version for an arbitrary volume. 
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Reynolds transport theorem for two arbitrary volumes A quantity Bu which is continuously 
defined over an arbitrary volume Vu is given by Bu = jy ^ dV. The volume Vu is assumed to 
consist of fictive mass and shall be material (conserving the fictive mass). Then the rate of 
change with time of Bu, according to the Reynolds transport theorem ([25]) . is given as 

^ [ ^pdV = [ ^dV + (f ifu-ndA , (26) 



dt / ^ J dt 

V, 



^ ^\pdV = [ ^dV + I ■dA+ I ipgu- {-n)dA . (28) 



dt y y 9t 

where u is the velocity of the boundary of the considered volume. 

Applying the general formula to the two control volumes together with (I26p yields Reynolds 
transport theorems for each of the two volumes Vi and Vg as 

^ [ ^dV = [ ^dV+ l\pv-dA+ [ ifiu-ndA, (27) 



and 



dij J dt 

Vg Vg Ag A 

By adding ()27p and (j28p we obtain a Reynolds transport theorem for the entire volume 
Vo = Vi + Vga.s 

^j\pdV = J ^dV + ^ ifv -dA- J^ipju-ndA . (29) 
Vo A 

Equation (|29p is the Reynolds transport theorem for a material volume with a singular in- 
terface. It states that the rate of accumulation of a quantity in a material volume, where ip 
undergoes a jump on an interface can be interpreted as the rate of accumulation of the quan- 
tity in a volume that equals the material volume at a given time plus convective flow of the 
quantities ipi and ipg leaving the volume through the outer surface and the interface at that 
time. Here again the integration limits are dropped where the integrals are evaluated at 
a given time. Only the integration domain of the integral along the interface A has to be 
indicated. 
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By using Gauss theoreir0, the Reynolds transport theorem (I29p for a material volume with 
a singular interface can be rewritten as 



J <fdV = j (^ + '^- b'^]) + y [^(i?--u)] -ndi. (30) 



Vo A 



3.3 Kinematics of a moving interface and velocities 

In this section we present the kinematics of the moving interface and discuss the material (or 
convected) surface derivative, fluid velocity and interface velocity. 

Kinematics of the moving interface If we consider phase change, the interface is not com- 
posed of a fixed set of particles, there will be mass transfer between the interface and the two 
adjoining phases. According to the basic assumption of continuum theory a surface consists 
at every moment of infinitely many particles. In particular, a particle joining the interface 
coincides with the particle that was at that position before. For a particle which is leaving the 
interface instantaneously another particle emerges. So although there is a many-to-one map- 
ping between particles and the region in the surface that is occupied by them, we can assign 
one representing particle for all possible particles at one point. We will call the representing 
particle simply particle (see [TrToj ) . 

The position vector to a point on the surface was given in section [2] as a function of surface 
coordinates and time ([T|) and is here denoted by a lower index S 

x^=Xg{u°',t) with a = 1,2. (31) 

At a given time a particle on the surface is represented by a position vector, which is here 
also denoted by a lower index S 

^s=^s(^o)- (32) 

We take this position as the initial position and call it intrinsic surface reference configuration. 
Conversely at a given time every position in the surface corresponds to a particle 

n^ = u^{U. (33) 



*For a material volume with an internal interface Gauss theorem becomes (compare footnote |3]) 
J'V-[(pv] dV = fipv-dA - Jjlipjv-ndA. 
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Obviously the initial position of a particle in the surface can be identified either by (1320 or 
by ()33p . At another time the particle is given by another set of coordinates as a function of 
the reference configuration of the particle and time 

u" = M"K,t). (34) 

With the assumption of a representing surface particle we established a one-to-one mapping 
between the coordinates of a surface particle and the surface coordinates, so that we can write 
reversely 

7x^ = n^(n",t) . (35) 

Equations ()34p and ()35p describe the intrinsic motion of a surface particle within the surface, 
without knowing how the surface itself is moving. 

The motion of a surface particle in space we get from the motion of the surface ()3ip and 
the intrinsic motion of the surface particles on the surface (I34p . (I35p as 

Xs=XsiuQ,t) or = Xg{u'^{u'^,t),t) . (36) 

Equation (j36P is not reversible in general. A position in space is corresponding to every surface 
particle, but the converse is not true (as explained for phase change). 

A surface field variable is here denoted by (pg. It can be given with pip as a function of 
space and time 

Or it can be given with ()36p as a function of particle and time 

Material surface derivative The partial derivative of a surface field variable with respect to 
time where the material surface coordinates are held constant is called material (or convected) 
surface derivative and it is here denoted with 



\ , lip 

It is the rate of change of a surface field variable with respect to time an observer moving 
with a surface particle would see. The material surface derivative of a position vector is the 
velocity of a given surface particle 

V = . (37) 
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Here we did not denote v with a lower index S to be consistent with the balance equations 
as they will be given later. Surface balance equations are often conveniently given in spatial 
surface coordinates. To obtain the surface material derivative of a surface field variable 
iPg{u°{u^,t),t) in spatial surface coordinates the chain rule has to be applied 

^sfs ^ 9yjg 9(/?g Dn° 
Bt dt 9n° Bt 

With the surface gradient pTj) 

and the intrinsic surface velocity 

w = a« (38) 
the material surface derivative becomes 

+ wV,.iPs. (39) 



Bt dt 

The material surface derivative is the local rate of change at a position of a given surface 
particle at a given time plus the convective rate of change related to the moving surface. 

Velocity of an interface particle relative to the velocity of a moving interface The material 
surface derivative of the surface position vector Xg{u°'{uQ ,t),t) is by (j39p given as 

^ = ^ + ^.V,.,. (40) 

The partial derivative of the surface position vector with respect to time (where u"' held 
constant) is the velocity of the moving interface 

u = ^. (41) 

For the second term on the right hand side we have [w ■ V) = w (V^ ® x^^ and the surface 
gradient (|lip of the surface position vector is the surface identity tensor 

- — ®x„ = or ® 
so that the material derivative of the surface position vector becomes 



V_g (g> Xg = a" - — (g) a^g = a" (g) = a'^ ® aa = I s : (42) 



= ^ + -VA. (43) 

V = u + w . 



17 



The intrinsic surface velocity is the velocity of a surface particle relative to the velocity of the 
surface 



w = V 



u 



(44) 



Note that in general u has a normal and a tangential part, so that w is not necessarily the 
tangential part of v. 



Surface velocity for an implicit surface parametrization In section 12.31 we discussed the 
geometrical properties of a surface defined by an implicit function F{Xg{u°^ ,t),t) = 0. Differ- 
entiating F = with respect to time gives 

dF 

'dt ~ 

where we dropped the subscript S at for simplicity. 



dF dF dx' 
'dt 9x* dt 



and equivalently 



u-VF = , 



(45) 



With the normal vector n = as derived in (jl5p we can write either 



u ■ n 



dF 
dt 



VF 



or 



u ■ n 



u 



, . (46) 
|VF| |VF| ^ ^ 

The first equation is independent of the parametrization, so that all possible surface velocities 
have the same normal component u ■ n, which is called the speed of displacement. It is 
convenient to choose a parametrization such that the surface velocity becomes the surface 
normal velocity 



u 



n (Xin It 



[u • n]n . 



The surface defined by F{u, v, z,t) = z — h{u, v, t) has the surface position vector 

Xg (n, v,t) = uSx + V ey + h(u, v, t) . 

For an implicit surface parametrization the surface velocity is the surface normal velocity and 
is given by 



u 



dx^ 
dt 






dh 
. dt. 



(47) 



Dotting (j43l) with n gives 



V. ■ n 



u ■ n -\ — ;r— otrw • n 



Dt 



u ■ n . 
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3.4 Reynolds transport theorem and divergence theorem for a surface 

In the balance equations we also need the Reynolds transport theorem for an interface, which 
is not material if we allow phase change. However, we can always assume the interface to be 
composed of a fixed set of particles, as in section [331 For a quantity Sq continuously defined 
over such an interface Aq, we write Sq = J^^ fgdA. The rate of change of Sq with respect to 
time is given by 
dSoit) d 



dt dt 

In an Eulerian representation AQ{Xg,t) depends on time. As for the Reynolds transport 
theorem for a material volume the area element is transformed with the surface Jacobian 

j from spatial coordinates into material coordinates dA = j dAo . 

Then Aq{^^) does not depend on time and integration and differentiation can be interchanged, 
so that the time derivative becomes the material surface derivative (1391) 



dtj I m 

Ao 



• 



Using the material surface derivative of the surface Jacobian = j Vg • we get 



After transforming the area element back into spatial coordinates and by using the material 
surface derivative ([5^]) with the relative velocity we get for the rate of change of Sq with 
respect to time 



dt 



+ [v,-u]-V,ip,+ip,V -v^] dA, (48) 



9^ 



This is the Reynolds transport theorem for surfaces. Comparing the last equation with the 
Reynolds transport theorem for material volumes ()24p shows an additional term related to 
the moving surface. 
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The last equation of (1480 can be transformed further. In this section we derive the di- 
vergence theorem for surfaces, which wih be used here and again in section 14.31 so that we 
derive it using the abbreviation / = ip^Vg. Sphtting / in normal and tangential part the 
area integral of the surface divergence of / becomes 



/di 



[{f -71)71] di + 



[(/•&) b] di. 



(49) 



Applying Product Rule to the first integrand on the right hand side gives 
V,-[{f -71) n] = n-V,{f ■ 71) + {f-7i)V,- 71 = -2Hf- 71, 



where we invoked the definition of mean curvature p3|) in section and the fact that the 
surface gradient Q is perpendicular to the normal vector on the surface. 

The second integral is an intrinsic area in- 
tegral. The binormal vector b = t x n is 
perpendicular to both the tangential vector 
t along the curve and to the surface normal 
vector, and therefore is the outward normal 
to the boundary curve as shown in figure HI 
Using Stokes theorenj§ to transform the area 
integral into a line integral Figure 4: Base vectors on a bounded surface 




equation 



s-[{f-m dA 

becomes 



{f-h)dC. 



V,-fdA 



2H f ■7idA+ i f ■ bdC . 



(50) 



This equation is called the divergence theorem for surfaces. (It should not be confused with 
the divergence theorem in footnote [5j) 



Alternative version of the Reynolds transport theorem for surfaces Substituting (f50|) into 
the last equation of (|48p gives an alternative version of the Reynolds transport theorem for 
surfaces 

^ J di = j (^^-u-V,^,-2H^^v^-7?j dA+j> ^^v^-bdC. (51) 



^Stokes theorem entirely defined in surface vectors is given as: j Vg ■ / AA = § f ■ bAC. 
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The rate of accumulation of a quantity in a surface can be interpreted as the rate of accumu- 
lation of the quantity in a material surface that equals the surface at a given time plus flux 
arising from the moving surface, plus convective flux normal to the surface (curvature term) 
and convective flux through the boundary curve of the area. 

4 Generic model equations for two-phase flows with surface tension 

In this section we derive the generic model equations for a problem with a moving interface. 
First we recall the balance equations for incompressible fluids and then we derive a differential 
balance equation for a moving interface between two fluids, a so called jump condition. This 
interface balance equation includes phase change, but to also include surface tension, a balance 
equation for the interface itself has to be formulated, which is then added to the balance 
equation for the bulk phases. 

The material about balance equations in general is mainly based on [BiStLi] and |S199] . 
Further on [De] and |Hu] . The material about interface balance equations is mainly based 
on |TrTo| , |S107| and [Edj . Some recent references on two phase flow problems are |IsHi] , |Ko| 
and |GaPr] . 

4.1 Balance equations for bulk fluids 

Balance equations are formulated for physical quantities that are continuously defined over a 
spatial region (for instances a volume), such as mass, momentum or energy. We denote those 
quantities by B = jyip dV. 

Balance equation for a material volume A material (mass conserving) volume Vq is in 
general moving with time. A balance equation for a physical quantity Bq = Jy^ ip dV in a 
material volume states that the rate of accumulation of the quantity in the material volume 
is given by an conductive flux of the quantity (not connected to mass) that enters the volume 
across the surface plus supply of the quantity to the material volumqj 



''Some authors distinguish between supply of quantity to the volume and production of quantity within the 
volume. Then conservation equations can be defined as balance equations without a production term. 
However, it is more intuitive to distinguish only between surface terms and volume terms. 




(52) 
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where C is the efflux density and vr is the supply density. The surface element vector dA = 
ndA is directed outwards, normal to the surface. As before we drop the integration limits 
except in the case, where the derivative of the integral is taken. With the Reynolds transport 
theorem (j24p the balance equation (j52p for a material volume becomes 

(53) 



Balance equation for a stationary volume To derive a balance equation for a stationary 
volume the Reynolds transport theorem in the form of ()25p is substituted for the right hand 
side of §21 

^dV + (f ifv - dA = -/ c-dA+ [ -KdV . 



ipdV = - 6 ipv ■ dA- 6 C- dA+ irdV . (54) 



9t „ 

Now the integration domain of the volume integral on the left hand side is constant and 
differentiation and integration can be interchanged. This gives the balance equation for a 
stationary volume 

d 

dt 

This equation has again a physical meaning: The accumulation of B = J ip dV in a stationary 
volume is given by convective and conductive flux of quantity across the surface to the volume 
plus supply of quantity to the volume. 

Differential balance equation For numerical computations a differential version of the bal- 
ance equation is preferable. The starting point is equation ()53p 

Using the divergence theorem (see footnote [3] in section 13. ip for the surface integral, we get 



/ (^ + '^ - + -vr)dy = 0. 



This equation must hold for any arbitrary volume so that we get the differential balance 
equation 

^ + V-{ipv) = -V-C + 7^ ■ (55) 

Although this equation represent the same physical phenomenon as before (accumulation, 
flux, supply) the various terms cannot be interpreted in the same way as the integral balance 
equations. 
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4.2 Jump conditions at an interface between two fluids 

Next we derive balance equations for an interface between two homogeneous bulk phases with 
phase change (condensation or evaporation). 

Balance equation for a material volume with a singular interface The balance equation 
(j52p for a quantity Bq in a material volume 



ifdV = - C-dA+ / TrdV 



_d 

dt 

Vo 

holds for a material volume whether or not the volume has a singular interface. The rate of 
accumulation in Vq = Vi + Vg is the sum of the rate of accumulation in the volumes Vi and Vg 

Vo Vl Vg 

The volumes Vi and Vg are not material, so we use Reynolds transport theorem for an arbitrary 
volume with a singular interface ([30]) . Substituting ([30|) for the left hand side of ([52]) gives 
the balance equation for a material volume with a singular interface 

/ (^^'^' ^'^'"^) ^^^/ [^^ ■'>^dA = - j) C-dA + J ndV . (56) 

A 

The double brackets denote again the difference between the two values at the surface. 



To derive a differential form of ([56]) . a spe- 
cial volume in a form of a small box is con- 
sidered, which is moving together with the 
interface as shown in figure [3 Two faces of 
the box are parallel to the interface. By tak- 
ing the limit of the shorter side faces Ah — t- 
the volume integrals vanish and Ai and Ag 
merge with A. For the volume integrals to 
vanish their integrands must be limited (but 
not necessary continuous). Then the nor- 
mal vectors of the two outer faces of the box 




Figure 5: Volume in form of a box 
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transform into either n or —n and only the surface integral over the interface A remains 



n + 



n 



dA = 0. 



The integral must hold for any arbitrary surface so that the integrand must be zero 

<f[v-u]-nj+ C-nj=0 . (57) 



Equation (l57l) is called a jump condition and describes the phase change at the interface and 
the conductive flux across the interface |S107] . 

This jump condition describes phase change and conductive flux across the surface, but it 
does not allow for modeling intrinsic surface properties such as surface tension (we took the 
limit of Ah^O). 



4.3 Balance equation and jump condition including surface tension 

Based on the assumption that a fluid interface is actually a three-dimensional region with a 
thickness of maybe one or more molecule diameters, the effect of the interface on the adjoining 
bulk fluids can be represented by assuming a two-dimensional interface consisting of surface 
excess mass. Surface mass is assumed to have similar properties as three-dimensional mass, 
such as surface density, surface viscosity, surface tension and so on. Then, analogous to the 
balance equation of three dimensional continua, a balance equation for the interface can be 
given. Adding this derived interface balance equation to the balance equation for a material 
volume with a singular interface gives a balance equation which includes surface tension. 



Balance equation for a surface The interface between condensate and vapor is not material, 
the fluid velocity differs from the velocity of the surface. Nevertheless a balance equation 
similar to the balance equation (I52p for a material volume can be given, as explained in 
section [3^ For a quantity Sq continuously defined over a surface A we write Sq = dA. 
Then a balance equation for Sq states that the rate of accumulation of surface quantity in 
the surface A is given by conductive flux of surface quantity across the boundary curve of the 
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surface plus supply of surface quantity at the surface 

^,dA = - j>Cs-'^C + j T,,dA, (58) 

A 

where C,^ is the surface flux density and tt^ is the surface supply density. The line element 
vector dC = m dC is directed outwards normal on the boundary curve, see figure H] in 
section [331 Using Reynolds theorem for a surface ([jS]) . equation ([58]) becomes 



+ • ^s] - ^ • "^s^s ] dA = - (bC,-dC+ TT^dA. (59) 



' " S LY-S ^SJ ' S^t^S I Y "*S ' I "S 

Differential balance equation for a surface To derive a differential version of (f59]) we trans- 
form the line integral into an area integral using the surface divergence theorem (|50p 

/ + ■ ""^^ " " ■ ^^'^^ + [V^ • + 2^Cs • '^l - ^s) = , 

where -ff is the mean curvature. This equation must hold for any arbitrary area so that the 
differential surface balance equation is given as 

^ + V, • [ifs Vs]-u. V,yp, = -[v^.C, + 2HCs-n]+7r^ . (60) 

Balance equation including phase change and surface tension Adding the surface balance 
equation ()59p to the balance equation for a material volume with a singular interface (|56p 
gives a balance equation that includes phase change and surface tension 



C-dA- Cs ■ dC + j TTdV + J TTgdA . (61) 



Jump condition including phase change and surface tension From (j6T]) we derive a jump 
condition in the same way as discussed in the last section. Transforming the line integral into 
an area integral using the surface divergence theorem (j50p and considering a small volume 
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enclosing the interface which we let shrink into a surface gives 



If [v — u]- n 



+ 



C-n 



dt 



If there is no material accumulation (no mass) in the surface, the surface density variables 
and VTg are zero and from the right hand side, only the underlined terms remain: 



(62) 



This jump condition describes conductive flux across the interface, phase change and surface 
tension. We made no additional assumption other than assume a continuous surface, in 
particular we do not allow the interface to break off. 

5 Model equations with phase change and surface tension 

In this section the balance equations for mass, momentum and energy for two incompressible 
fluids and for the interface between them are obtained and we discuss appropriate simplifica- 
tions. In the last section we summarize the deduced system of partial differential equations 
including the jump conditions and discuss boundary conditions for the system of partial dif- 
ferential equations. 



5.1 Mass, momentum and energy equation 

Mass The mass balance equation for an incompressible fluid is given by (j55p with ip = p 
and C = vr = 



^ + V.(pt;) = 0, 
which gives with the assumption of constant density 
V • = . 



(63) 



(64) 



If the Mach number of a fluid is small compared to unity the fluid can be considered as an 
incompressible fluid [BiStLij . 

For momentum and energy another balance equation which makes use of (j63p is more 
preferable. For that we substitute <p in (j55]) by pil), apply product rule on both terms on the 
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left side and receive 
9(pV) 



dt 



+ V-[{p -0) v] 



dp 

dt 



+ V -(pv) 



+P 



"9^ 



+ V ■ Vip 



The first bracket is zero according to (j63p . The second bracket is the material derivative of 
as derived in section 13.11 

Thus we get the generic balance equation ()55p in an equivalent form 

Dip 



P 



Bt 



Momentum The momentum equation we get by substituting in (j65p ip 
TT = pg (where ip and tt are vectors and ^ is a second order tensor) 



P 



9? + ^ 



V -S + pg 



(65) 



-S and 



(66) 



The stress tensor can be divided into a contribution of the fluid at rest and the fluid in motion, 
S = —pI + T. The body force vector is the gravity vector g assuming that there are no other 
body forces. 

If the fluids are Newtonian (linearly viscous) fluids the (deviatoric part of the) stress tensor 
for both phases is given by T = [Vi» + (V?;)"^] + y (V • v), with shear viscosity p and 
modified bulk viscosity p' . Together with the incompressibility condition the momentum 
equation 



P 



9? + ^ 



— Vp + pV V + p g 



then forms the well known Navier-Stokes equations. 

Note that if we had derived the momentum equation from (j55p , the divergence term in the 
convective term would be nonlinear. The advantage of deriving the momentum equation from 
(j65p is that the divergence term is then linear which makes numerical discretization easier. 
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Energy According to the first law of thermodynamics the increase of internal and kinetic 
energy in a material control volume is given by heat supply plus power due to work acting on 
the fluid. The differential equation governing the energy is 



D 



(V • [Sv] + pg-v) + (- V ■q + pz) 



(67) 



where e is the internal energy per unit mass, Sv is the power due to surface forces per unit 
area, pg v is the gravity power per unit volume, q the heat flux per unit area and z the heat 

2 

production per unit volume, which is zero in our case. Here ■0 = e + C, = —Sv + q and 
vr = pg-v + p z. To get the energy in a more commonly used form, we subtract the mechanical 
energy equation from (|67|) . The mechanical energy equation is formed by the scalar product 
of momentum equation and velocity. By this we ge10 



P 



9e 



Ve 



-'V ■ q + S 'Vv 



Assuming that the heating effect of friction can be neglected, the dissipative term T • Vt? is 
zero. Moreover, S • 'Vv vanishes since S • 'Vv = T • 'Vv — p(I • 'Vv) = —p (V ■ v) = 0. 

Constitutive equations for internal energy and heat flux complete the equations. For small 
temperature differences internal energy can be described by a linear function e = c(T — Tq) + 
e(To), where c is the specific heat capacity. The heat flux is given by Fourier's law q = — AVT, 
where A is the heat conductivity. This results in the heat equation 
^9T 



pc 



dt 



+ V-VT 



AV^T 



The material properties viscosity, heat capacity and heat conductivity are in general func- 
tions of density, pressure and temperature, but for incompressible fluids only temperature 
dependency need to be considered. If the temperature interval between the two phases is 
small the material properties can be assumed to be constant. 



5.2 Mass, momentum and energy jump conditions 

Mass We start with (j62p . Assuming the bulk variables '■p = p and C = 0, and the surface 
variable Co = 0) gives for the mass jump condition at the interface 



p [v — u]- n 







Equation ()68p states that the amount of mass flux that enters the interface must also leave 
the interface m/ = nig (so we can omit the index). 



''Here we used the identity V ■ {Sv) ^ {V ■ S) ■ v + S • Vv 
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Momentum To get the momentum jump condition at the surface between condensate and 
vapor, we set in (f62]l = pv, C, = —S and C,^ = —S^. As for the stress tensor, the surface 
stress tensor can be divided into two components S ^ = al ^ + T ^ , where I g is the surface 
identity tensor ([9]). Assuming a surface without mass we have Tg = and from the surface 
stress tensor, only the interfacial tension remains Sg = alg, where a is the surface tension 
coefficient, see section 14.31 Surface tension can be understood as the counterpart of the 
pressure in the bulk fluid. With this substitution the momentum jump condition becomes 



pv [v — u] ■ n 



Sn 



W + 2H an . 



(69) 



If we neglect temperature dependency of the surface tension coefficient (no Marangony effects) 
(J is constant within the surface. 

We split the vector equation (j69p into three scalar equations by multiplying it first with the 
normal vector and then with the two tangential vectors. The tangential equations are equal, 
so we skip the third equation and use the symbol t to denote both tangential vectors. Further 
we make use of (j68p and the assumption of no-slip at the surface \ v • t =0. Then we get 



rhv ■ n + p — n • Tn 



t ■ Tn 



2Ha 







(70) 
(71) 



aI,+T,)u 



2 

Energy The energy jump condition is with = p (e+^), C = —Sv+q, = ( 
—au and by making the same assumptions as for the momentum jump condition given by 

„.2 



P[e + 



V — u 



n 



n ■ Sv 



+ 



q ■ n 



((ju) + 2 H a n ■ u 



(72) 



where e is the internal energy. Instead of a balance equation for the internal energy a for- 
mulation with the enthalpy h = e + ^ is more convenient, because enthalpy is a measurable 
quantity whereas internal energy is not easy to measure. 

Writing n ■ Sv = v ■ Sn (S is symmetric), splitting the vectors on the left side of (|72p into 
their normal and tangential components according to a - b = (a ■ n){b ■ n) + (a ■ t){b ■ t) , and 
applying the chain rule to Vg • [a uj gives 



(v ■ n)2 (v ■ tf 



{v ■n){n- Sn) + {v ■ t) {t ■ Sn) 



+ 



qn 



a'Vg-u + u-'Vc;Cr + 2Han-u. (73) 
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Subtracting from the energy jump condition ()73p the scalar product of surface velocity u = 
(u ■ n)n and momentum jump condition ()69p 

m {u ■ n){v ■ n) — (u ■ n) (n • Sn) | = w • [Vg a + 2H a ■ n] 

gives together with the no-slip condition and the tangential momentum jump condition ([7T 
after rearranging the pressure term 



P 

m \ e-\ — 
P 



(V ■ n] 



{u ■ n){v ■ n)^ — — ti] • n (n • '^i^)'^ + |9 " ''^| 



a Vg • u . 



Using again Chain Rule and that n • Vg = we get for the right hand side 

(T Vg • u = o" Vg ■ {u ■ n) n = a {u ■ n) Vg ■ n = 2a H u ■ n 
Expanding the kinetic energy term 



r {v ■ nf - 




I 2 J 





{v ■ n — u ■ n) 



+ {v ■ n){u ■ n) 



{u • n)^ 



and noting that the surface velocity jump is zero gives the energy jump condition as 
{[v — u] ■ n)^ 



mh\\ + 



m ■ 



\v — u] ■ n (n ■ Tn) 



+ 



q-n = 2a H u ■ n . 



(74) 



If kinetic energy, viscous energy and the effect of surface tension can be neglected, the energy 
jump condition for an interface between two fluids becomes 

m A/i+ q • n| = , (75) 
where A/i = /i| is the latent heat of vaporization. 
5.3 Summary of the model equations 

The two fluids (i = l,g) are described by the continuity equation for an incompressible fluid, 
the momentum equation and the energy equation: 



Continuity equation 



const, density 
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Momentum equation 

inertia 



Energy equation 

ft Ci ( — + • VTi 

transient + convection 



pressure friction gravity 



heat conduction 



Newtonian fluid 
const, viscosity 



no dissipation 
const, heat conductivity 
const, thermal capacity 



The jump conditions at tlic interface between the two fluids for mass, momentum and 
energy are with |(^| := ipg — ipi (dropping the tilde on the normal and tangential vectors): 



Mass 

|p [t) — u] • n 
Momentum 



mass flux 



m v • nj + |pj — |n ■ T • n 



momentum 
due to 

condensation 



surface pressure 



const, surface tension coefficient 
no slip between both phases 

surface tension 



t • T • n] 

shear stress 



= 



Energy 



no kinetic energy 
no dissipation 



m Ah 

condensation 



q ■ n 

heat flux 



with T = n [Vv + (Vv)^] and q = -AVT. 
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Boundary conditions As many boundary conditions for each coordinate of an unknown are 
necessary as the equation has derivatives of this unknown. One more equation is missing and 
for this we take the condition of thermodynamic equilibrium at the interface, according to 
which the fluid temperatures are equal at the moving surface 

Ti = Tg. (76) 

Prom the five jump conditions we need one equation to calculate the mass flux, so that 
four equations remain to calculate the boundary conditions for three velocity components, 
pressure and temperature. We can use either the mass jump condition or the energy jump 
condition to calculate the mass flux. 

From the remaining equations either the mass jump condition or the normal momentum 
jump condition can be used to calculate one velocity boundary condition, depending on which 
equation is used to compute the mass flux. Which of the two equations is used depends on 
further simplifications of the problem. Often the term with the mass flux in the normal 
momentum jump condition is small and can be dropped, then this equation is not available 
to calculate the mass flux. 

The second and third velocity boundary condition are taken from the tangential momentum 
jump conditions. 

SI Units of some variables and material properties 
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H 


m ^ 


mean curvature 


S 


N 


stress tensor (J = kgms~^) 


T 


N 


viscous stress tensor 


T 


K 


temperature (C = K - 273,15) 


c 


Jkg-iR-i 


specific heat capacity ( J = N m) 


e 


Jkg-i 


specific internal energy 


9 


_2 

ms 


gravity vector {g = 9,81 ms"^) 


a 


Wm-^K"^ 


heat transfer coefficient (W = Js~^) 


Ah 


Jkg-i 


latent heat of vaporization 


rh 


kgs~^ 


volume specific mass fiux 


P 




pressure 


Q 




heat flux vector 


u,v 




velocity 


1^ 


kg m"^ s"^ 


dynamical viscosity 


A 


Wm-^K-i 


thermal conductivity 


P 


kg 


density 


a 


N 


surface tension 
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